Exploration of real bulk RNA-seq datasets

TB datasets

Introduction

This is a curated dataset containing the Africa, Brazil, India, and the G6 datasets, which are all TB datasets with active TB patients, patients that will progress to active TB, and people with TB infection who don’t progress. There is a batch effect especially with the Brazil data. This object is a SummarizedExperiment with 7 assays. The ‘combat’ assay was made by running the logtpm values through combat with ‘SequencingBatch’ as the batch and ‘Label’ as the covariate. “Label” is the condition Non-progressor, Progressor, and Active.

Overall design

The combined dataset contains 4 studies, in which Brazil dataset has 2 batches within the study:

• Africa: 197 samples, (104 Non-progressor, 77 Progressor, 16 Active)
• Brazil: 51 samples, (21 Non-progressor, 16 Progressor, 14 Active)
  • Batch 1: 39 samples, (9 Non-progressor, 16 Progressor, 14 Active)
  • Batch 2: 12 samples, (12 Non-progressor, 0 Progressor, 0 Active)
• G6: 399 samples, (304 Non-progressor, 95 Progressor, 0 Active)
• India: 44 samples, (19 Non-progressor, 0 Progressor, 25 Active)

Note: Gender is an important factor in this study.

Load data & summarize basic info

After loading the data, we calculate a few basic statistics of the dataset.

######  Samples  #######
- Total number of samples:
[1] 691

- Number of samples in condition groups:
group
Non-progressor     Progressor         Active 
           448            188             55 

- Number of samples in batches:
batch
  Africa Brazil_1 Brazil_2       G6    India 
     197       39       12      399       44 

In each batch separately:
         Non-progressor Progressor Active
Africa              104         77     16
Brazil_1              9         16     14
Brazil_2             12          0      0
G6                  304         95      0
India                19          0     25


######  Genes  #######
- Total number of genes:
[1] 24391

- Number of genes with 0 counts:
14504 59.5%

- Number of genes with maximum 10 counts:
2587 10.6%

- Number of genes with ONLY 0s:
0 0%


######  Counts  #######
- Range of count in whole dataset:
[1]       0 7674206

- Average count in whole dataset:
[1] "Mean: 2350.9"
[1] "Median: 75"

There are no genes with only zero values across all samples (probably because of cleaning in the curation stage). The figure below shows distribution of percentages of zero counts for the remaining genes (each data point represents a gene, i.e., how many zeros are there in each gene).

######  Proportion of genes with X% zeros across samples  #######
 <10.0%  <20.0%  <30.0%  <40.0%  <50.0%  <60.0%  <70.0%  <80.0%  <90.0% <100.0% 
  0.646   0.033   0.023   0.020   0.020   0.021   0.025   0.029   0.045   0.139 

Existence of batch effect

In this merged dataset, there are both “study effects” across the 4 studies (targeting different populations), and batch effect within the same study (Brazil dataset).

PCA

PCA (Brazil only)

Dendrograms (Brazil only)

Library size

Next, we take a look at the observed library sizes of samples, and to see if there is a difference in library size across batches.

Range of observed library sizes across samples:
[1] "4.91e+05" "1.09e+08"

Gene expression distribution

Our next step is to interrogate the gene expression distributions.

Fold changes in mean

We use edgeR to perform differential expression analysis, and calculate (log-) fold changes of all genes between condition groups and batches.

Change WRT condition

To rule out the impact of batch effect, we limit our DE analysis with regard to biological condition in one of the batches. Here we specify the batch to be G6 due to a larger batch size, and evaluate DE between progressors and non-progressors.

####  Differential genes relative to biological condition 
####  (within batch 1 / batch specified by 'test_batch')
test_batch <- "G6"
counts_batch <- assay(rds_obj,"counts")[, batch==test_batch]
cond_batch <- as.factor(as.character(group[batch==test_batch]))
y <- DGEList(counts=counts_batch, group=as.factor(cond_batch))
y <- calcNormFactors(y, method="TMM")
design <- model.matrix(~as.factor(cond_batch))
y <- estimateDisp(y, design)

######  Changes in terms of average gene counts  ######
               top50 up top50 down 100 least DE
Non-progressor   1728.5     3362.2        757.1
Progressor       3117.4     2773.3        706.5

######  Range of fold change across biological conditions  ######
          top50 up top50 down 100 least DE
min FC       1.501      0.004        0.934
median FC    3.975      0.705        1.000
mean FC     21.299      0.563        0.998
max FC     668.720      0.899        1.000

######  Range of Log(fold change) across biological conditions  ######
             top50 up top50 down 100 least DE
min logFC       0.406     -5.553       -0.068
median logFC    1.380     -0.350        0.000
mean logFC      3.059     -0.574       -0.002
max logFC       6.505     -0.107        0.000

Conclusions:

• Biological fold change comparing average counts in top DE genes is about 1.2-1.8.
• Gene-wise fold change in top DE genes has a wide range (median 4, mean 21).

Change WRT batch

We are also interested to see how the average (mean) gene counts change across batches (i.e. degree of mean batch effect).

To explore this mean batch effect, we first perform DE with regard to the batch variable. To rule out the impact of biological condition, we perform this analysis within Non-progressors only, and compare Africa & G6 studies. Note that this DE is not true DE - it’s testing if there is a significant mean difference between two batches.

####  Differential genes relative to batch 
####  (within condition 0 / condition specified by 'test_cond')
rds_subset <- rds_obj[, rds_obj$Dataset %in% c("Africa", "G6")]
test_cond <- "Non-progressor"
counts_cond <- assay(rds_subset, "counts")[, rds_subset$Label==test_cond]
batch_cond <- as.factor(as.character(rds_subset$SequencingBatch[rds_subset$Label==test_cond]))
y <- DGEList(counts=counts_cond, group=as.factor(batch_cond))
y <- calcNormFactors(y, method="TMM")
design <- model.matrix(~as.factor(batch_cond))
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)
qlf <- glmQLFTest(fit, coef=2)
de_res <- topTags(qlf, n=nrow(rds_subset))
up_genes <- rownames(de_res)[de_res$table$logFC >= 0]  # all up-regulated genes
head_up_genes <- head(up_genes, n=50)  # top-50 up-regulated genes
tail_up_genes <- tail(up_genes, n=50)  # 50 least up-regulated genes
batch_up_inds <- up_genes
down_genes <- rownames(de_res)[de_res$table$logFC < 0]  # all down-regulated genes
head_down_genes <- head(down_genes, n=50)  # top-50 down-regulated genes
tail_down_genes <- tail(down_genes, n=50)  # 50 least down-regulated genes
batch_down_inds <- down_genes

######  Batch differences in terms of mean gene counts  ######
       top50 up top50 down 100 least changed
Africa  1917.33    4100.75             96.39
G6      6824.55    3466.82            142.55

######  Range of fold change between batches  ######
- All genes - FC:
0.005 1291.191

- All genes - log(FC):
-5.352 7.163

- Top genes - FC:
          top50 up top50 down 100 least DE
min FC       1.872      0.005        0.953
median FC    4.424      0.425        1.000
mean FC     52.612      0.410        0.996
max FC    1291.191      0.609        1.000

- Top genes - log(FC):
             top50 up top50 down 100 least DE
min logFC       0.627     -5.352       -0.049
median logFC    1.487     -0.856        0.000
mean logFC      3.963     -0.892       -0.004
max logFC       7.163     -0.496        0.000

Conclusions:

• Mean batch effect fold change comparing average counts in top genes ranges from 1.18-3.56. Remember that these are genes most significantly different across batch.
• Gene-wise mean batch effect fold change in top genes has a wide range (mean 4.4, median 52.6).

Batch mean

Median (of all genes) gene-wise mean:

               Africa     G6 Brazil1
Non-progressor  87.85 129.77   15.00
Progressor      74.19 119.88   22.75

Scaled by that of Africa:

               Africa   G6 Brazil1
Non-progressor      1 1.48    0.17
Progressor          1 1.62    0.31

Batch dispersion

Aside from average gene counts, we are interested in comparing dispersions in the batches. Difference in the dispersions across batches suggests variance batch effect not accounted for by Poisson model. We visualize the gene-wise dispersion estimates from edgeR:

Design matrix not provided. Switch to the classic mode.
####  Statistics about dispersion estimates from the two batches  ####
               Africa       G6   India Brazil_1 Brazil_2
min disp.      0.0054   0.0034  0.0001   0.0390   0.0305
median disp.   0.1916   0.2422  0.1616   0.6837   0.2086
mean disp.     3.3203   2.7797  0.3211   3.0653   0.4466
max disp.    102.4000 102.4000 50.4797  38.5716   6.2874

Conclusions:

• Mean dispersion in the Africa / Batch 1 in the Brazil study is about 10 times of those in the India study.
• We see huge max gene-wise dispersions, possibly caused by 0 counts.

Batch variance

Median (of all genes) gene-wise variance:

                Africa      G6 Brazil1
Non-progressor 1967.72 4722.60  135.19
Progressor     1620.59 4623.17  469.53

Scaled by that of Africa:

               Africa   G6 Brazil1
Non-progressor      1 2.40    0.07
Progressor          1 2.85    0.29

Impact of 0 counts on gene distributions

We have seen above that over half of genes in the dataset contain some portions of zeros. Having 0 counts can have a good impact on mean and dispersion estimation. We evaluate such impact in this section, and see how NB distribution fits to non-zero portion of genes. We first take the non-zero portion of the genes.

Fold change in mean in non-zero portion

Change WRT condition in G6

Within G6, we compare mean (here refers to non-zero mean) across biological condition groups. Specifically, we saved the previous DE analysis results using edgeR, on counts containing zeros. We call all genes with \(\log(FC) >= 0\) as “up” genes, and \(\log(FC) < 0\) as “down” genes. The top 50 “up” and “down” genes are still determined by the previous DE containing zeros.

######  Changes in MEAN non-zero gene counts  ######
               top50 up top50 down 100 least DE
Non-progressor   2169.9     4497.1       2196.1
Progressor       3914.2     3786.5       2085.1

######  Changes in MEDIAN non-zero gene counts  ######
               top50 up top50 down 100 least DE
Non-progressor      609       2236        474.5
Progressor          961       1834        471.0

Conclusions:

• Biological fold change (in non-zero portion) comparing average counts in top genes is about 1.2-1.8 (mean).
• O counts did not strongly affect biological signal in terms of average counts in top genes .

Change WRT batch in non-progressors

Similarly, we compare non-zero mean in top genes across batches within “Non-progressors” of Africa and G6 studies. Here “top genes” refer to top gene lists from previous DE within control samples comparing batches.

######  Changes in MEAN non-zero gene counts  ######
       top50 up top50 down 100 least DE
Africa   2072.8     4110.2        446.7
G6       6834.0     3530.7        609.6

######  Changes in MEDIAN non-zero gene counts  ######
       top50 up top50 down 100 least DE
Africa      473       2413           30
G6         1600       2054           34

Conclusions:

• Mean batch effect fold change (non-zero portion) comparing average counts in top genes ranges from 1.16-3.3 (mean).
• O counts did not strongly affect mean batch effect in terms of average counts in top genes.

Batch mean with 0 counts removed

Median (of all genes) gene-wise mean after removing zeros:

               Africa    G6 Brazil1
Non-progressor 118.62 155.7   35.33
Progressor     109.31 173.2   33.94

Scaled by that of Africa:

               Africa   G6 Brazil1
Non-progressor      1 1.31    0.30
Progressor          1 1.58    0.31

Batch dispersion with 0 counts removed

Dispersion estimates are more likely impacted by 0 counts. We re-estimate the dispersions in each batch, using non-zero portion of each gene. Then we compare the dispersion estimates across batches.

####  Statistics about dispersion estimates (from NON-ZERO portions) from the two batches  ####
             Africa     G6  India Brazil_1 Brazil_2
min disp.    0.0000 0.0000 0.0000   0.0000   0.0000
median disp. 0.1293 0.1739 0.1176   0.5544   0.1708
mean disp.   0.1885 0.2359 0.1601   0.5460   0.2548
max disp.    5.0019 4.1344 3.1316   3.7155   4.8316

Conclusions:

• Mean dispersion of Brazil batch 1 is roughly 3.4 times that of India.
• The absolute value of max gene-wise dispersion are reduced after 0 counts are removed.
• In this report, dispersion are calculated by different methods before and after removing 0 counts. Not directly comparable.
• We show in another analysis (TB_run_nonzero_disps.R) that the difference in dispersion is reduced after removing 0s.

Batch variance with 0 counts removed

Median gene-wise variance after removing zeros:

                Africa       G6 Brazil1
Non-progressor 4871.43  8110.23  964.11
Progressor     4561.47 11948.42 1571.00

Scaled by that of Africa:

               Africa   G6 Brazil1
Non-progressor      1 1.66    0.20
Progressor          1 2.62    0.34

Is there a difference in zero portions across batches?

Another question related to batch effect, is whether there is any difference in the percentage of 0 counts in a gene between two batches.

Conclusions:

• There isn’t an obvious difference in zero fractions between batches (aside from Brazil batch 2 which is too small).

Goodness of fit for NB distribution to non-zero portion of genes

It is difficult to evaluate quantitatively the goodness-of-fit of NB distribution on data. Instead, we generate Q-Q plots for 100 randomly selected genes, comparing data with NB distribution using ML estimated parameters. We generate these plots within G6/Africa, for cases and control samples separately. Figures are stored under {output_dir}/qqplots_{Study}/.

More justification for having different dispersion parameter for different batches

Finally, a very important assumption of ComBat-seq model (compared to adding batch covariates in linear model of edgeR and DESeq), is that having a different dispersion parameter for each batch is beneficial for DE performance. As argued by many, if not all, of these papers, accurate estimation of dispersion paramters are important for DE detection. As we have seen above, there is some difference in dispersion distribution across batches. If dispersion estimates within each batch and using the whole study are not consistent with each other, then it is justification that we may need different dipsersion parameter for different batches, and it could be good for DE.

Here we generate scatter plots comparing dispersion estimates within each batch to estimates from whole dataset (including batch as a covariate in whole study).

In Brazil

Gene-wise dispersion estimates

Trended dispersion estimates

Africa VS G6

Gene-wise dispersion estimates

Trended dispersion estimates

Brazil Batch 1 VS India

Gene-wise dispersion estimates

Trended dispersion estimates

Session info

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] MASS_7.3-51                 ggdendro_0.1-20             scales_1.0.0                gridExtra_2.3              
 [5] reshape2_1.4.3              ggplot2_3.0.0               dendextend_1.9.0            edgeR_3.22.5               
 [9] limma_3.36.5                DESeq2_1.20.0               SummarizedExperiment_1.10.1 DelayedArray_0.6.6         
[13] BiocParallel_1.14.2         matrixStats_0.54.0          Biobase_2.40.0              GenomicRanges_1.32.7       
[17] GenomeInfoDb_1.16.0         IRanges_2.14.12             S4Vectors_0.18.3            BiocGenerics_0.26.0        

loaded via a namespace (and not attached):
 [1] bitops_1.0-6           bit64_0.9-7            RColorBrewer_1.1-2     rprojroot_1.3-2        prabclus_2.2-6        
 [6] tools_3.5.1            backports_1.1.2        R6_2.3.0               rpart_4.1-13           Hmisc_4.1-1           
[11] DBI_1.0.0              lazyeval_0.2.1         colorspace_1.3-2       trimcluster_0.1-2.1    nnet_7.3-12           
[16] withr_2.1.2            tidyselect_0.2.5       bit_1.1-14             compiler_3.5.1         htmlTable_1.12        
[21] labeling_0.3           diptest_0.75-7         checkmate_1.8.5        DEoptimR_1.0-8         mvtnorm_1.0-8         
[26] robustbase_0.93-3      genefilter_1.62.0      stringr_1.3.1          digest_0.6.18          foreign_0.8-71        
[31] rmarkdown_1.10         XVector_0.20.0         base64enc_0.1-3        pkgconfig_2.0.2        htmltools_0.3.6       
[36] htmlwidgets_1.3        rlang_0.3.0            rstudioapi_0.8         RSQLite_2.1.1          bindr_0.1.1           
[41] jsonlite_1.5           mclust_5.4.1           acepack_1.4.1          dplyr_0.7.7            RCurl_1.95-4.11       
[46] magrittr_1.5           modeltools_0.2-22      GenomeInfoDbData_1.1.0 Formula_1.2-3          Matrix_1.2-14         
[51] Rcpp_0.12.19           munsell_0.5.0          viridis_0.5.1          stringi_1.2.4          whisker_0.3-2         
[56] yaml_2.2.0             zlibbioc_1.26.0        flexmix_2.3-14         plyr_1.8.4             grid_3.5.1            
[61] blob_1.1.1             crayon_1.3.4           lattice_0.20-35        splines_3.5.1          annotate_1.58.0       
[66] locfit_1.5-9.1         knitr_1.20             pillar_1.3.0           fpc_2.1-11.1           geneplotter_1.58.0    
[71] XML_3.98-1.16          glue_1.3.0             evaluate_0.12          latticeExtra_0.6-28    data.table_1.11.8     
[76] gtable_0.2.0           purrr_0.2.5            kernlab_0.9-27         assertthat_0.2.0       xtable_1.8-3          
[81] viridisLite_0.3.0      class_7.3-14           survival_2.42-6        tibble_1.4.2           AnnotationDbi_1.42.1  
[86] memoise_1.1.0          bindrcpp_0.2.2         cluster_2.0.7-1       

References